###checkpoint_cor
###gene与免疫检查点相关性分析#气泡图
#引用包
setwd("E:\\1.Ovarian_Treg\\1.picture\\9.checkpoint")
library(limma)
library(reshape2)
library(ggplot2)
library(ggpubr)
expFile="OV_TPM.txt"      #表达输入文件
riskFile="rs_tcga.txt"       #风险输入文件
geneFile="gene.txt"       #免疫检查点的基因文件

#读取基因表达文件,并对数据进行处理
rt=read.table(expFile, header=T, sep="\t", check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)

#读取基因文件
gene=read.table(geneFile, header=F, sep="\t", check.names=F)
sameGene=intersect(row.names(data),as.vector(gene[,1]))
data=t(data[sameGene,])
data=log2(data+1)

#删除正常样品
group=sapply(strsplit(row.names(data),"\\-"),"[",4)
group=sapply(strsplit(group,""),"[",1)
group=gsub("2","1",group)
data=data[group==0,]
row.names(data)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*?)\\-.*","\\1\\-\\2\\-\\3",row.names(data))
data=avereps(data)

#合并数据
risk=read.table(riskFile, sep="\t", header=T, check.names=F, row.names=1)
sameSample=intersect(row.names(data),row.names(risk))
rt1=cbind(data[sameSample,],risk[sameSample,])
rt1=rt1[,c(sameGene,"riskScore")]

data=rt1

#读取风险文件
sameSample=intersect(row.names(data), row.names(risk))
data=data[sameSample,,drop=F]
data=data[,-ncol(data)]

model_gene=read.table("modelgene.txt", sep="\t", header=T, check.names=F)
data_model=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data_model=avereps(data_model)
data_model=data_model[model_gene[,1],]
data_model=t(data_model)

group=sapply(strsplit(row.names(data_model),"\\-"),"[",4)
group=sapply(strsplit(group,""),"[",1)
group=gsub("2","1",group)
data_model=data_model[group==0,]
row.names(data_model)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*?)\\-.*","\\1\\-\\2\\-\\3",row.names(data_model))
data_model=avereps(data_model)
sameSample=intersect(row.names(risk), row.names(data_model))
data_model=data_model[sameSample,,drop=F]
risk=risk[sameSample,,drop=F]
data_model_risk=cbind(data_model,risk[,"riskScore",drop=F])
data_model_risk=data_model_risk[row.names(data),]
risk=data_model_risk
#相关性分析
outTab=data.frame()
for(checkpiont in colnames(data)){
  for(gene in colnames(risk)){
    x=as.numeric(data[,checkpiont])
    y=as.numeric(risk[,gene])
    corT=cor.test(x,y,method="spearman")
    cor=corT$estimate
    pvalue=corT$p.value
    text=ifelse(pvalue<0.001,"***",ifelse(pvalue<0.01,"**",ifelse(pvalue<0.05,"*","")))
    outTab=rbind(outTab,cbind(Gene=gene, checkpiont=checkpiont, cor, text, pvalue))
  }
}
outTab=outTab[outTab$text%in%c("***","**","*"),]
save(outTab,file = "cor_checkpoint.rdata")
# 添加正负相关信息：
data_result=outTab
Neg_Pos <- rep("Neg",184)
Neg_Pos[which(data_result$cor>0)] <- "Pos"

data_result$Neg_Pos <- Neg_Pos
#data_result$pvalue <- c(rep(0.00002,5),rep(0.0002,5),runif(180,0,0.08))

# 修改p值表示方法：
pvalue <- rep(NA,184)
data_result$pvalue=as.numeric(data_result$pvalue)
pvalue[which(data_result$pvalue > 0.05)] <- ">0.05"
pvalue[which(data_result$pvalue < 0.05)] <- "<0.05"
pvalue[which(data_result$pvalue < 0.01)] <- "<0.01"
pvalue[which(data_result$pvalue < 0.001)] <- "<0.001"
pvalue[which(data_result$pvalue < 0.0001)] <- "<0.0001"
table(pvalue)
data_result$pvalue <- pvalue
unique(data_result$Gene)
data_result$cor=as.numeric(data_result$cor)
# 绘图：
data_result$Gene=factor(data_result$Gene,levels = unique(data_result$Gene))

# 绘图：
ggplot(data_result,aes(checkpiont,Gene))+
  # 蓝色气泡图：
  geom_point(aes(fill=pvalue, size=abs(cor)),color = "#999999",shape=21)+
  scale_fill_manual(values = c("#212c5f","#3366b1","#42b0e4","#7bc6ed","#dfe1e0"))+
  # 红色气泡图：
  geom_point(data = data_result[which(data_result$Neg_Pos == "Pos"),],
             aes(color=pvalue, size=abs(cor)),shape=16)+
  scale_color_manual(values = c("#f26666","#f49699","#facccc","#facccc","#d9dbd9"))+
  # 主题：
  theme_bw()+
  theme(panel.grid.minor.x = element_blank(),
        panel.grid.major.x = element_blank(),
        # 坐标轴label方向：
        axis.text.x = element_text(angle = 45, hjust = 1),
        # 图例间距：
        legend.margin = margin(20,unit = 'pt'),
        legend.title = element_text(size = 8),
        legend.text = element_text(size = 8))+
  xlab("")+
  ylab("")+
  # 图例：
  guides(size = guide_legend(title = "Spearman's p"),
         fill = guide_legend(title = expression("Negtive \ncorrelation \np-value")),
         col = guide_legend(title = expression("Positive \ncorrelation \np-value")))

ggsave("cor_buble_heatmap1.pdf",height = 3,width = 12)
